{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Application in which we read and process altimetry measurements using `diva`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import re\n",
    "import glob\n",
    "import logging\n",
    "import numpy as np\n",
    "from matplotlib import rcParams\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import pydiva2d\n",
    "import netCDF4\n",
    "import divaaltimetry\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a `ModuleNotFoundError` message, ensure that the module `pydiva2d` is available in the Python path.\n",
    "\n",
    "# Setup \n",
    "\n",
    "## Logging configuration\n",
    "\n",
    "The *logging* is already configured in `pydiva2d`.<br>\n",
    "Replace 'DEBUG' by 'INFO', 'WARNING' or 'ERROR'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "logger = logging.getLogger('divaAltimetry')\n",
    "logger.setLevel(logging.DEBUG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Files and directories\n",
    "\n",
    "Set the path to the Diva installation you want and the main input files.   \n",
    "A different file will be used for the mesh generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "divadir = \"/home/ctroupin/Software/DIVA/DIVA-diva-4.7.1\"\n",
    "projectdir = \"/home/ctroupin/Projects/Altimetry-Interpolation/\"\n",
    "coastfile = os.path.join(projectdir, \"diva/coast.cont\")\n",
    "datadir = os.path.join(projectdir, \"data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/\")\n",
    "paramfile = os.path.join(projectdir, \"diva/param.par\")\n",
    "paramfilemesh = os.path.join(projectdir, \"diva/param.par.mesh\")\n",
    "outputdir = os.path.join(projectdir, \"results/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva2D/\")\n",
    "figdir = \"../figures/diva2d/\"\n",
    "if not os.path.exists(divadir):\n",
    "    logger.error(\"Diva directory doesn't exist\")\n",
    "if not os.path.exists(outputdir):\n",
    "    os.makedirs(outputdir)\n",
    "    logger.debug(\"Create output directory\")\n",
    "if not os.path.exists(figdir):\n",
    "    os.makedirs(figdir)\n",
    "    logger.debug(\"Create figure directory\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matplotlib options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rcParams.update({'font.size': 14, 'figure.dpi': 300, 'savefig.bbox': 'tight'})\n",
    "coordinates = (-6.75, 36.001, 30, 48.)\n",
    "meridians = np.arange(-8., 40., 8.)\n",
    "parallels = np.arange(30., 50., 4.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Projection\n",
    "Will be used for the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = Basemap(projection='merc',\n",
    "            llcrnrlon=coordinates[0], llcrnrlat=coordinates[2],\n",
    "            urcrnrlon=coordinates[1], urcrnrlat=coordinates[3],\n",
    "            lat_ts=0.5 * (coordinates[2] + coordinates[3]), resolution='i')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Directories\n",
    "Create variables storing the Diva directories and files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DivaDirs = pydiva2d.DivaDirectories(divadir)\n",
    "DivaFiles = pydiva2d.Diva2Dfiles(DivaDirs.diva2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis\n",
    "## Create mesh\n",
    "As the domain is always the same (surface level), we only create the mesh once at the beginning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of nodes: 6704\n",
      "Number of interfaces: 16293\n",
      "Number of elements: 9496\n"
     ]
    }
   ],
   "source": [
    "mesh2d = pydiva2d.Diva2DMesh().make(divadir, contourfile=coastfile, paramfile=paramfilemesh)\n",
    "mesh2d.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diva calculation\n",
    "We use another parameter file, loop over the data files and copy the outputs (netCDF) in the specified directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140201_20140221.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140202_20140222.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140203_20140223.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140204_20140224.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140205_20140225.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140206_20140226.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140207_20140227.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140208_20140228.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140209_20140301.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140210_20140302.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140211_20140303.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140212_20140304.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140213_20140305.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140214_20140306.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140215_20140307.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140216_20140308.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140217_20140309.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140218_20140310.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140219_20140311.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140220_20140312.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140221_20140313.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140222_20140314.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140223_20140315.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140224_20140316.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140225_20140317.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140226_20140318.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140227_20140319.dat\n",
      "/home/ctroupin/Projects/Altimetry-Interpolation/data/SEALEVEL_MED_PHY_L3_REP_OBSERVATIONS_008_049/diva/data_20140228_20140320.dat\n"
     ]
    }
   ],
   "source": [
    "datafilelist = sorted(glob.glob(os.path.join(datadir, 'data_201402*.dat')))\n",
    "logger.info(\"Working on {} files\".format(len(datafilelist)))\n",
    "\n",
    "datepattern = re.compile(r'data_(\\d{8})_(\\d{8})')\n",
    "\n",
    "for datafile in datafilelist:\n",
    "    \n",
    "    print(datafile)\n",
    "    outputfile = \"\".join((os.path.basename(datafile).split('.')[0], '.nc'))\n",
    "    figname = os.path.basename(datafile).split('.')[0]\n",
    "    \n",
    "    match = datepattern.search(os.path.basename(datafile))\n",
    "    if match:\n",
    "        d1 = \"-\".join((match.group(1)[:4], match.group(1)[4:6], match.group(1)[6:8]))\n",
    "        d2 = \"-\".join((match.group(2)[:4], match.group(2)[4:6], match.group(2)[6:8]))\n",
    "        figtitle = \"$-$\".join((d1, d2)) \n",
    "    logger.info(\"Output file: {0}\".format(outputfile))\n",
    "    \n",
    "    results2d = pydiva2d.Diva2DResults().make(divadir, datafile=datafile,\n",
    "                                              paramfile=paramfile, \n",
    "                                              contourfile=coastfile,\n",
    "                                              outputfile=os.path.join(outputdir, outputfile))\n",
    "    \n",
    "    \n",
    "    # Make plot\n",
    "    SLA = divaaltimetry.AltimetryField().from_diva2d_file(os.path.join(outputdir, outputfile))\n",
    "    SLA.add_to_plot(figname=os.path.join(figdir, figname), figtitle=figtitle, m=m,\n",
    "                    meridians=meridians, parallels=parallels,\n",
    "                    vmin=-0.2, vmax=0.2,\n",
    "                    cmap=plt.cm.RdYlBu_r)    \n",
    "    \n",
    "logger.info(\"Summary:\")\n",
    "logger.info(\"Figures available in {}\".format(os.path.abspath(figdir)))\n",
    "logger.info(\"Results available in {}\".format(os.path.abspath(outputdir)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Diva-python3.6",
   "language": "python",
   "name": "diva-python3.6"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
